marshall_visualizations
Emma Rieves & Zac Clement
Created 4/27/2022, Updated 7/5/2022
Load the Libraries
library(sf)## Warning: package 'sf' was built under R version 4.1.3## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUElibrary(sp)## Warning: package 'sp' was built under R version 4.1.3library(gstat)## Warning: package 'gstat' was built under R version 4.1.3library(raster)## Warning: package 'raster' was built under R version 4.1.3library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.1.3## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --## v ggplot2 3.3.6 v purrr 0.3.4
## v tibble 3.1.7 v dplyr 1.0.9
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.1.3## Warning: package 'tibble' was built under R version 4.1.3## Warning: package 'tidyr' was built under R version 4.1.3## Warning: package 'readr' was built under R version 4.1.3## Warning: package 'dplyr' was built under R version 4.1.3## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()library(spacetime)## Warning: package 'spacetime' was built under R version 4.1.3library(stringr)
library(readr)
library(leaflet)
library(gridExtra)##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combinelibrary(gifski)## Warning: package 'gifski' was built under R version 4.1.3library(nngeo)## Warning: package 'nngeo' was built under R version 4.1.3# to make the spherical geometry errors go away
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched offImport Data
Air quality & sensor data
# AIR QUALITY
## corrected air quality
AQ_df = read.csv("../intermediary_outputs/corrected_AQ_data.csv")
## sensor info
sensors = read.csv("../intermediary_outputs/sensor_data_full.csv")Geometries
# BOUNDARIES
# Boulder boundaries
BO_CO = st_read("../GIS_inputs_destruction_fireboundary/Boulder_county_munis/Municipalities.shp")## Reading layer `Municipalities' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Boulder_county_munis\Municipalities.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -105.525 ymin: 39.91405 xmax: -105.0528 ymax: 40.23771
## Geodetic CRS: WGS 84fire_counties = BO_CO %>%
dplyr::filter(ZONEDESC == "Louisville" |
ZONEDESC == "Superior" |
ZONEDESC == "Broomfield" |
ZONEDESC == "Lafayette" |
ZONEDESC == "Boulder")
# Boulder Boundary
boulder_precincts = st_read("../GIS_inputs_destruction_fireboundary/Unincorporated_Boulder/Unincorporated_Boulder.shp")## Reading layer `Unincorporated_Boulder' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Unincorporated_Boulder\Unincorporated_Boulder.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 129 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -105.4297 ymin: 39.91371 xmax: -105.0528 ymax: 40.21256
## Geodetic CRS: WGS 84# Broomfield Boundary
broomfield_precincts = st_read("../GIS_inputs_destruction_fireboundary/Broomfield_Precincts/Precincts.shp")## Reading layer `Precincts' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Broomfield_Precincts\Precincts.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 61 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 3093759 ymin: 1202697 xmax: 3150881 ymax: 1259386
## Projected CRS: NAD83(HARN) / Colorado North (ftUS)# Westminster Boundary
westminster_city = st_read("../GIS_inputs_destruction_fireboundary/Westminster_CityLimits/CityLimits.shp")## Reading layer `CityLimits' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\Westminster_CityLimits\CityLimits.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 18 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -105.1654 ymin: 39.81811 xmax: -104.9877 ymax: 39.96854
## Geodetic CRS: WGS 84# make all the precincts have the same columns (we don't need a lot of the data)
boulder_precincts <- boulder_precincts %>%
select(c(OBJECTID, geometry, SHAPEarea, SHAPElen))
broomfield_precincts <- broomfield_precincts %>%
mutate(SHAPEarea = Shape__Are,
SHAPElen = Shape__Len,
OBJECTID = OBJECTID + 1000) %>% # add 1,000 to the broomfield precincts to get rid of overlapping labels
select(-c("GlobalID", "GIS_ID", "PRECINCT_N", "MAP_COLOR", "USER_COMME", "LAST_UPDAT", "Shape__Are", "Shape__Len"))
westminster_city <- westminster_city %>%
mutate(SHAPEarea = ShapeSTAre,
SHAPElen = ShapeSTLen) %>%
select(c(OBJECTID, geometry, SHAPEarea, SHAPElen))
# save proj string for fire-affected counties to use for other transformations
prg = raster::crs(fire_counties,asText=TRUE)
# all boundaries joined together
all_bounds <- boulder_precincts %>%
rbind(st_transform(broomfield_precincts, crs=prg)) %>%
rbind(st_transform(westminster_city, crs=prg))
# remove the holes that existed so we just have the outside boundary of this area
all_bounds <- all_bounds %>%
st_union() %>%
st_remove_holes()## although coordinates are longitude/latitude, st_union assumes that they are planar# Marshall fire boundary
wfigs_fire = st_read("../GIS_inputs_destruction_fireboundary/WFIGS_-_Wildland_Fire_Perimeters_Full_History/FH_Perimeter.shp")## Reading layer `FH_Perimeter' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\WFIGS_-_Wildland_Fire_Perimeters_Full_History\FH_Perimeter.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 9500 features and 108 fields (with 4 geometries empty)
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -165.8152 ymin: 13.38021 xmax: 144.6906 ymax: 69.05174
## Geodetic CRS: WGS 84# filter fire to marshall fire only; update crs
marshall_fire = wfigs_fire %>% filter(poly_Incid == "Marshall") %>% st_transform(crs = st_crs(prg))After getting the boundary data, we need to read in the geographic locations of buildings that were damaged or destroyed by the fire.
# BUILDINGS
# Read in destroyed home and building sites
destroyed_homes = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/destroyed_homes.shp")## Reading layer `destroyed_homes' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\destroyed_homes.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1030 features and 25 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -105.23 ymin: 39.93075 xmax: -105.135 ymax: 40.01238
## Geodetic CRS: WGS 84damaged_homes = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/damaged_homes.shp")## Reading layer `damaged_homes' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\damaged_homes.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 143 features and 25 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -105.2309 ymin: 39.93081 xmax: -105.1439 ymax: 40.01083
## Geodetic CRS: WGS 84destroyed_businesses = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/destroyed_businesses.shp")## Reading layer `destroyed_businesses' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\destroyed_businesses.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 24 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -105.1651 ymin: 39.94129 xmax: -105.1384 ymax: 39.95788
## Geodetic CRS: WGS 84damaged_businesses = st_read("../GIS_inputs_destruction_fireboundary/output_damage_files/damaged_businesses.shp")## Reading layer `damaged_businesses' from data source
## `C:\Users\zclem\Documents\GitHub\marshall_fires\GIS_inputs_destruction_fireboundary\output_damage_files\damaged_businesses.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 24 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -105.2132 ymin: 39.94129 xmax: -105.1384 ymax: 40.01153
## Geodetic CRS: WGS 84# create 'type' column; select only the columns needed; make projection string match
destroyed_homes = destroyed_homes %>%
mutate(type = "destroyed - residential") %>%
dplyr::rename(jurisdiction = JURISDICTI) %>%
dplyr::select(jurisdiction, type, latlong, geometry) %>%
st_transform(crs = st_crs(prg))
destroyed_businesses = destroyed_businesses %>%
mutate(type = "destroyed - non-residential") %>%
dplyr::rename(jurisdiction = Jurisdicti) %>%
dplyr::select(jurisdiction, type, latlong, geometry) %>%
st_transform(crs = st_crs(prg))
damaged_homes = damaged_homes %>%
mutate(type = "damaged - residential") %>%
dplyr::rename(jurisdiction = JURISDICTI) %>%
dplyr::select(jurisdiction, type, latlong, geometry) %>%
st_transform(crs = st_crs(prg))
damaged_business = damaged_businesses %>%
mutate(type = "damaged - non-residential") %>%
dplyr::rename(jurisdiction = Jurisdicti) %>%
dplyr::select(jurisdiction, type, latlong, geometry) %>%
st_transform(crs = st_crs(prg))
# combine into one df, using the "," indicator keeps decimal points
destroyed_damaged = rbind(destroyed_homes,
destroyed_businesses,
damaged_homes,
damaged_business) %>%
separate(latlong, c("lat","long"), ",", convert = FALSE)
# change lat and long to floats
destroyed_damaged$lat = as.double(destroyed_damaged$lat)
destroyed_damaged$long = as.double(destroyed_damaged$long)Leaflet Maps
# create "location status" column for mapping]
sensors$loc_status = str_c(sensors$Location, " - ", sensors$Status)
# make sure we only have valid sensors that have data
sensors = na.omit(sensors)# set factors to false
options(stringsAsFactors = FALSE)
# create palette for destroyed and damaged building colors for leaflet map
pal = colorFactor(c("red","orange","magenta","light pink"), domain = unique(destroyed_damaged$type))
## create indoor/outdoor icons
sensorIcon = awesomeIconList(
"inside" = makeAwesomeIcon(
icon = "home",
iconColor = "white",
markerColor = "blue",
library = "fa"
),
"outside" = makeAwesomeIcon(
icon = "tree",
iconColor = "white",
markerColor = "green",
library = "fa"
))
sensors## X ID fire_period post_fire_period pre_fire_period
## 1 1 3062 100.000000 100.00000 100.00000
## 2 2 4583 100.000000 99.93082 100.00000
## 3 3 4643 100.000000 99.33126 100.00000
## 4 4 15663 100.000000 99.52727 100.00000
## 5 5 27120 48.333333 99.86740 100.00000
## 6 6 30743 74.000000 73.96770 100.00000
## 7 7 33109 53.333333 99.70022 100.00000
## 8 8 38623 7.333333 95.18621 100.00000
## 10 10 45325 99.666667 88.22783 100.00000
## 11 11 45369 99.666667 99.94235 100.00000
## 13 13 50133 100.000000 100.00000 100.00000
## 14 14 53729 100.000000 99.84434 100.00000
## 16 16 56209 100.000000 94.93255 96.66667
## 17 17 56431 86.666667 99.87893 100.00000
## 18 18 60507 100.000000 99.93188 100.00000
## 19 19 60517 100.000000 99.82705 100.00000
## 20 20 64671 100.000000 99.81552 100.00000
## 21 21 64967 22.333333 99.86740 100.00000
## 22 22 65669 17.666667 99.89046 100.00000
## 23 23 66173 100.000000 99.96541 100.00000
## 24 24 67917 99.333333 99.96541 100.00000
## 25 25 68375 7.333333 99.25055 100.00000
## 26 26 69103 4.000000 99.72328 100.00000
## 28 28 70241 100.000000 99.96541 100.00000
## 29 29 70911 100.000000 94.94408 100.00000
## 31 31 72931 100.000000 99.94811 100.00000
## 34 34 79121 7.333333 92.70725 100.00000
## 35 35 79907 38.000000 99.99423 100.00000
## 36 36 81149 100.000000 99.59068 100.00000
## 37 37 82879 100.000000 99.68292 100.00000
## 38 38 83423 100.000000 98.81817 100.00000
## 39 39 83597 97.333333 99.58492 100.00000
## 40 40 83661 100.000000 95.07091 100.00000
## 41 41 83785 4.000000 99.72328 100.00000
## 42 42 84169 100.000000 99.91929 100.00000
## 43 43 84585 76.000000 98.76780 100.00000
## 44 44 84697 1.333333 96.62170 100.00000
## 45 45 84777 10.666667 63.97674 100.00000
## 46 46 85643 100.000000 100.00000 100.00000
## 47 47 85999 14.000000 97.89000 100.00000
## 48 48 86477 13.666667 98.16672 100.00000
## 49 49 86637 99.666667 99.78669 100.00000
## 50 50 87465 100.000000 99.98847 100.00000
## 51 51 87531 100.000000 99.38314 100.00000
## 52 52 88035 100.000000 99.52150 100.00000
## 53 53 88533 100.000000 99.97694 100.00000
## 55 55 90255 3.666667 99.64257 100.00000
## 56 56 90297 89.000000 90.30900 100.00000
## 57 57 91877 100.000000 99.94811 100.00000
## 58 58 92321 100.000000 99.97117 100.00000
## 59 59 94411 100.000000 88.86199 100.00000
## 60 60 95591 100.000000 100.00000 100.00000
## 61 61 95831 99.000000 87.35155 100.00000
## 62 62 96731 100.000000 100.00000 100.00000
## 63 63 97647 100.000000 99.93658 100.00000
## 64 64 100521 100.000000 99.60798 100.00000
## 65 65 101003 4.000000 99.40044 100.00000
## 66 66 102614 100.000000 99.47538 100.00000
## 67 67 102724 100.000000 99.82128 100.00000
## 69 69 103724 7.333333 98.20708 100.00000
## 70 70 103820 49.000000 46.41992 100.00000
## 71 71 104594 99.666667 99.78093 100.00000
## 72 72 106214 22.333333 99.85587 100.00000
## 73 73 106492 48.666667 98.85276 100.00000
## 74 74 106800 100.000000 99.67139 100.00000
## 75 75 106906 99.666667 99.99423 100.00000
## 76 76 107466 63.000000 99.96541 100.00000
## 78 78 108558 99.666667 99.77516 100.00000
## 80 80 110042 100.000000 100.00000 100.00000
## 81 81 110122 20.666667 99.96541 100.00000
## 82 82 111834 56.666667 99.34855 100.00000
## 83 83 112082 56.666667 99.59068 100.00000
## 84 84 112522 100.000000 97.51528 100.00000
## 85 85 112606 99.666667 97.58446 100.00000
## 86 86 112974 100.000000 99.98847 100.00000
## 87 87 113054 100.000000 99.96541 100.00000
## 88 88 113266 100.000000 99.87317 100.00000
## 93 93 115303 99.666667 99.94811 100.00000
## 94 94 117647 100.000000 99.82705 100.00000
## 95 95 119547 99.666667 99.44656 100.00000
## 97 97 120477 100.000000 99.96541 100.00000
## 99 99 120847 53.333333 99.98812 100.00000
## 101 101 121197 100.000000 99.98656 100.00000
## 103 103 122071 98.333333 99.00842 100.00000
## 106 106 122585 4.000000 96.54675 100.00000
## 115 115 124319 53.333333 99.46962 100.00000
## 200 200 141956 45.000000 64.57397 100.00000
## Status Month
## 1 Complete data throughout fire period Before or during 12-2021
## 2 Complete data throughout fire period Before or during 12-2021
## 3 Complete data throughout fire period Before or during 12-2021
## 4 Complete data throughout fire period Before or during 12-2021
## 5 Sensor offline during fire, returned online Before or during 12-2021
## 6 Sensor offline during fire, did not return online Before or during 12-2021
## 7 Sensor offline during fire, returned online Before or during 12-2021
## 8 Sensor offline during fire, returned online Before or during 12-2021
## 10 Complete data throughout fire period Before or during 12-2021
## 11 Complete data throughout fire period Before or during 12-2021
## 13 Complete data throughout fire period Before or during 12-2021
## 14 Complete data throughout fire period Before or during 12-2021
## 16 Complete data throughout fire period Before or during 12-2021
## 17 Complete data throughout fire period Before or during 12-2021
## 18 Complete data throughout fire period Before or during 12-2021
## 19 Complete data throughout fire period Before or during 12-2021
## 20 Complete data throughout fire period Before or during 12-2021
## 21 Sensor offline during fire, returned online Before or during 12-2021
## 22 Sensor offline during fire, returned online Before or during 12-2021
## 23 Complete data throughout fire period Before or during 12-2021
## 24 Complete data throughout fire period Before or during 12-2021
## 25 Sensor offline during fire, returned online Before or during 12-2021
## 26 Sensor offline during fire, returned online Before or during 12-2021
## 28 Complete data throughout fire period Before or during 12-2021
## 29 Complete data throughout fire period Before or during 12-2021
## 31 Complete data throughout fire period Before or during 12-2021
## 34 Sensor offline during fire, returned online Before or during 12-2021
## 35 Sensor offline during fire, returned online Before or during 12-2021
## 36 Complete data throughout fire period Before or during 12-2021
## 37 Complete data throughout fire period Before or during 12-2021
## 38 Complete data throughout fire period Before or during 12-2021
## 39 Complete data throughout fire period Before or during 12-2021
## 40 Complete data throughout fire period Before or during 12-2021
## 41 Sensor offline during fire, returned online Before or during 12-2021
## 42 Complete data throughout fire period Before or during 12-2021
## 43 Complete data throughout fire period Before or during 12-2021
## 44 Sensor offline during fire, returned online Before or during 12-2021
## 45 Sensor offline during fire, did not return online Before or during 12-2021
## 46 Complete data throughout fire period Before or during 12-2021
## 47 Sensor offline during fire, returned online Before or during 12-2021
## 48 Sensor offline during fire, returned online Before or during 12-2021
## 49 Complete data throughout fire period Before or during 12-2021
## 50 Complete data throughout fire period Before or during 12-2021
## 51 Complete data throughout fire period Before or during 12-2021
## 52 Complete data throughout fire period Before or during 12-2021
## 53 Complete data throughout fire period Before or during 12-2021
## 55 Sensor offline during fire, returned online Before or during 12-2021
## 56 Complete data throughout fire period Before or during 12-2021
## 57 Complete data throughout fire period Before or during 12-2021
## 58 Complete data throughout fire period Before or during 12-2021
## 59 Complete data throughout fire period Before or during 12-2021
## 60 Complete data throughout fire period Before or during 12-2021
## 61 Complete data throughout fire period Before or during 12-2021
## 62 Complete data throughout fire period Before or during 12-2021
## 63 Complete data throughout fire period Before or during 12-2021
## 64 Complete data throughout fire period Before or during 12-2021
## 65 Sensor offline during fire, returned online Before or during 12-2021
## 66 Complete data throughout fire period Before or during 12-2021
## 67 Complete data throughout fire period Before or during 12-2021
## 69 Sensor offline during fire, returned online Before or during 12-2021
## 70 Sensor offline during fire, did not return online Before or during 12-2021
## 71 Complete data throughout fire period Before or during 12-2021
## 72 Sensor offline during fire, returned online Before or during 12-2021
## 73 Sensor offline during fire, returned online Before or during 12-2021
## 74 Complete data throughout fire period Before or during 12-2021
## 75 Complete data throughout fire period Before or during 12-2021
## 76 Sensor offline during fire, returned online Before or during 12-2021
## 78 Complete data throughout fire period Before or during 12-2021
## 80 Complete data throughout fire period Before or during 12-2021
## 81 Sensor offline during fire, returned online Before or during 12-2021
## 82 Sensor offline during fire, returned online Before or during 12-2021
## 83 Sensor offline during fire, returned online Before or during 12-2021
## 84 Complete data throughout fire period Before or during 12-2021
## 85 Complete data throughout fire period Before or during 12-2021
## 86 Complete data throughout fire period Before or during 12-2021
## 87 Complete data throughout fire period Before or during 12-2021
## 88 Complete data throughout fire period Before or during 12-2021
## 93 Complete data throughout fire period Before or during 12-2021
## 94 Complete data throughout fire period Before or during 12-2021
## 95 Complete data throughout fire period Before or during 12-2021
## 97 Complete data throughout fire period Before or during 12-2021
## 99 Sensor offline during fire, returned online Before or during 12-2021
## 101 Complete data throughout fire period Before or during 12-2021
## 103 Complete data throughout fire period Before or during 12-2021
## 106 Sensor offline during fire, returned online Before or during 12-2021
## 115 Sensor offline during fire, returned online Before or during 12-2021
## 200 Sensor offline during fire, did not return online Before or during 12-2021
## Lon Lat Name Location
## 1 -105.2713 40.04077 Kalmia outside
## 2 -105.1680 40.13818 Clover Basin Consulting outside
## 3 -105.2515 39.98943 JENSA outside
## 4 -105.1550 40.10230 FRCC Niwot location outside
## 5 -105.2600 39.98123 Table Mesa Court outside
## 6 -105.1283 40.13157 Creekside outside
## 7 -105.2622 39.98305 Bear outside
## 8 -105.1749 39.98661 76th & S Bldr Rd outside
## 10 -105.1647 40.06149 Kincross Way outside
## 11 -105.1649 40.06169 Kincross Way inside
## 13 -105.1727 40.10403 HRNS-TEST inside
## 14 -105.0457 39.99003 Anthem Ranch outside
## 16 -105.2691 39.97224 Stony Hill Road outside
## 17 -105.2821 40.03435 Grape_and_Broadway outside
## 18 -105.2694 40.01557 Bella outside
## 19 -105.1199 39.87919 STANDLEY LAKE outside
## 20 -105.2709 40.04019 K2055i inside
## 21 -105.2458 39.98447 Majestic Heights, Boulder, CO outside
## 22 -105.1477 39.97686 Louisville outside
## 23 -105.1214 39.94361 Terracina Broomfield Apartments outside
## 24 -105.2833 40.00194 Chautauqua Heights outside
## 25 -105.1500 39.98119 Louisville Air Sensor outside
## 26 -105.3219 40.04352 Hawk Lane outside
## 28 -105.2804 40.02177 MKHOME2 outside
## 29 -105.2508 39.97368 Shanahan Ridge inside
## 31 -105.2279 39.99787 AZTEC CENTRAL outside
## 34 -105.1741 39.98653 Kitchen inside
## 35 -105.2716 40.03448 Vista Drive inside
## 36 -105.1549 39.97095 ntsky outside
## 37 -105.2401 39.98775 Mango Manor Coop inside
## 38 -105.2849 40.05431 Shining Mountain Waldorf School outside
## 39 -105.2563 39.97676 South Boulder outside
## 40 -105.1233 40.13444 Creekside and Sunset outside
## 41 -105.3218 40.04347 Hawk Lane inside
## 42 -105.2312 40.01628 Snowcave inside
## 43 -105.2308 39.99910 645 Manhattan Place inside
## 44 -105.2798 39.93268 Eldorado Springs outside
## 45 -105.1660 39.95349 Superior Town Hall outside
## 46 -105.3186 40.02700 Lower Seven Hills outside
## 47 -105.1554 39.93897 Superior - North Pool outside
## 48 -105.1552 39.92155 Superior - South Pool outside
## 49 -105.2515 39.97368 Viele Lake outside
## 50 -105.2845 40.06458 Dakota Ridge outside
## 51 -105.2394 39.98681 The Haunted Pickle Mansion inside
## 52 -105.1928 40.07116 Gunbarrel Green Spotted Horse (inside) inside
## 53 -105.1489 39.99006 Purple house indoors inside
## 55 -105.2551 39.98033 Miami Way outside
## 56 -105.0926 39.96745 South Pointe outside
## 57 -105.1451 39.98836 Eisenhower Dr outside
## 58 -105.1491 39.99013 Purple House outside
## 59 -105.2689 40.03923 Parkside outside
## 60 -105.2893 40.05224 Wonderland Lake outside
## 61 -105.1641 40.06173 Gunbarrel Hill outside
## 62 -105.2875 40.02247 Concord Ave outside
## 63 -105.0223 39.92575 Broomfield - Calkins Pl outside
## 64 -105.1932 40.07457 PurpleHome inside
## 65 -105.3151 40.04732 House inside
## 66 -105.0764 40.04520 Ponderosa outside
## 67 -105.1251 40.02127 Blue Heron Estates outside
## 69 -105.1499 39.97538 Home inside
## 70 -105.2813 40.00443 Casper’s House inside
## 71 -105.2007 40.07604 Bean Mountain Ln outside
## 72 -105.2458 39.98451 Majestic Heights inside
## 73 -105.3364 40.07924 Boulder Heights inside
## 74 -105.1249 40.02118 Blue Heron Estates inside
## 75 -105.2670 39.97523 1333 Wildwood Ct., Boulder, CO 80305 Inside inside
## 76 -105.2559 39.98311 Yale Road Home inside
## 78 -105.1720 40.09753 CountryCreek2 outside
## 80 -105.2562 40.03099 Valmont and 29th outside
## 81 -105.2427 39.98229 Hanover Ave. inside
## 82 -105.2451 40.09622 Snead Ct outside
## 83 -105.2549 40.09183 Eagle Ct. outside
## 84 -105.0474 39.82646 Linda Park Addition outside
## 85 -105.2670 39.97515 1333 Wildwood Ct., Boulder, CO 80305 outside
## 86 -105.1615 40.13421 Summerlin_Lane outside
## 87 -105.2837 40.04738 Quince Ave inside
## 88 -105.2803 40.04895 Redwood & Riverside outside
## 93 -105.0864 39.90192 Amli Arista outside
## 94 -105.1197 39.87945 Stanley Lake inside
## 95 -105.0439 39.95181 Broadlands outside
## 97 -105.1246 39.99547 Bohannan house outside
## 99 -105.2659 39.97925 Yard outside
## 101 -105.2376 40.00849 Hancock Drive outside
## 103 -105.1962 40.04964 EMSonJay outside
## 106 -105.2556 39.98821 Bates inside
## 115 -105.2567 39.98218 Hartford Drive inside
## 200 -105.1550 39.97581 Nighthawk Resident inside
## address
## 1 2065, Kalmia Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 2 5017, William Place, Longmont, Boulder County, Colorado, 80503, United States
## 3 3565, Eastman Avenue, Boulder, Boulder County, Colorado, 80305, United States
## 4 7037, Longview Drive, Boulder, Niwot, Boulder County, Colorado, 80503, United States
## 5 2799, Table Mesa Court, Boulder, Boulder County, Colorado, 80305, United States
## 6 2443, Eagleview Circle, Longmont, Boulder County, Colorado, 80504, United States
## 7 732, Ithaca Drive, Boulder, Boulder County, Colorado, 80305, United States
## 8 7618, South Boulder Road, Boulder, Boulder County, Colorado, 80303, United States
## 10 8059, Kincross Way, Heatherwood, Boulder, Boulder County, Colorado, 80301, United States
## 11 8059, Kincross Way, Heatherwood, Boulder, Boulder County, Colorado, 80301, United States
## 13 My Mama's Pie, Murray Street, Boulder, Niwot, Boulder County, Colorado, 80544, United States
## 14 4512, Silver Mountain Loop, Anthem, Broomfield, Colorado, 80032, United States
## 16 1916, Stony Hill Road, Boulder, Boulder County, Colorado, 80305, United States
## 17 Broadway & Grape Ave, Broadway, Washington Village, Boulder, Boulder County, Colorado, 80306, United States
## 18 1933, Grove Street, Boulder, Boulder County, Colorado, 80302, United States
## 19 10098, Oak Court, Walnut Grove, Westminster, Jefferson County, Colorado, 80021, United States
## 20 2099, Kalmia Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 21 605, South 42nd Street, Boulder, Boulder County, Colorado, 80305, United States
## 22 661, Cleveland Avenue, Louisville, Boulder County, Colorado, 80027, United States
## 23 Building 30, 13630, Via Varra, Overlook District, Broomfield, Colorado, 80020, United States
## 24 847, Cascade Avenue, Central Boulder - University Hill, Boulder, Boulder County, Colorado, 80302, United States
## 25 599, West Sagebrush Court, Louisville, Boulder County, Colorado, 80027, United States
## 26 67, Hawk Lane, Boulder County, Colorado, 80304, United States
## 28 2357, 13th Street, Washington Village, Boulder, Boulder County, Colorado, 80304, United States
## 29 Boulder, Boulder County, Colorado, 80305, United States
## 31 556, Aztec Drive, Boulder, Boulder County, Colorado, 80303, United States
## 34 South Boulder Road, Boulder, Boulder County, Colorado, 80027, United States
## 35 1970, Vista Drive, Boulder, Boulder County, Colorado, 80304, United States
## 36 775, West Lois Court, Louisville, Boulder County, Colorado, 80027, United States
## 37 4662, Ingram Court, Boulder, Boulder County, Colorado, 80305, United States
## 38 Shining Mountain Waldorf School, 999, Violet Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 39 1264, Ithaca Drive, Boulder, Boulder County, Colorado, 80305, United States
## 40 1839, Caleta Trail, Longmont, Boulder County, Colorado, 80504, United States
## 41 67, Hawk Lane, Boulder County, Colorado, 80304, United States
## 42 1707, Commerce Street, Boulder, Boulder County, Colorado, 80301, United States
## 43 645, Manhattan Drive, Boulder, Boulder County, Colorado, 80303, United States
## 44 Eldorado Springs Pool, Artesian Drive, Eldorado Springs, Boulder County, Colorado, 80025, United States
## 45 Superior Municipal Court, East William Street, Original Superior, Superior, Boulder County, Colorado, 80027, United States
## 46 442, Seven Hills Drive, Seven Hills, Boulder County, Colorado, 80304, United States
## 47 Town of Superior North Pool Facility, South Indiana Street, Rock Creek Ranch, Boulder, Superior, Boulder County, Colorado, 80027, United States
## 48 Superior South Pool, Huron Peak Drive, Rock Creek Ranch II, Boulder, Boulder County, Colorado, 80027, United States
## 49 3270, Everett Drive, Boulder, Boulder County, Colorado, 80305, United States
## 50 4949, Broadway, Boulder County, Colorado, 80304, United States
## 51 4682, Ingram Court, Boulder, Boulder County, Colorado, 80305, United States
## 52 5310, Spotted Horse Trail, Boulder, Boulder County, Colorado, 80301, United States
## 53 1998, Eisenhower Drive, Louisville, Boulder County, Colorado, 80027, United States
## 55 950, Miami Way, Boulder, Boulder County, Colorado, 80305, United States
## 56 South Point Drive, South Pointe, Lafayette, Boulder County, Colorado, 80026-2872, United States
## 57 370, Eisenhower Drive, Louisville, Boulder County, Colorado, 80027, United States
## 58 1998, Eisenhower Drive, Louisville, Boulder County, Colorado, 80027, United States
## 59 22nd Street, Boulder, Boulder County, Colorado, 80304, United States
## 60 501, Utica Avenue, Boulder, Boulder County, Colorado, 80303, United States
## 61 8067, Kincross Way, Heatherwood, Boulder, Boulder County, Colorado, 80301, United States
## 62 Concord Alley, Washington Village, Boulder, Boulder County, Colorado, 80302, United States
## 63 2830, Calkins Place, Broomfield, Colorado, 80020, United States
## 64 5583, Homestead Way, Boulder, Boulder County, Colorado, 80301, United States
## 65 230, Timber Lane, Pine Brook Hill, Boulder County, Colorado, 80304, United States
## 66 11821, Flatiron Drive, Leyner, Erie, Boulder County, Colorado, 80026, United States
## 67 2432, Ginny Way, Blue Heron Estates, Lafayette, Boulder County, Colorado, 80026, United States
## 69 681, Spruce Circle, Louisville, Boulder County, Colorado, 80027, United States
## 70 941, Lincoln Place, Central Boulder - University Hill, Boulder, Boulder County, Colorado, 80302, United States
## 71 6655, Bean Mountain Lane, Boulder, Boulder County, Colorado, 80301, United States
## 72 605, South 42nd Street, Boulder, Boulder County, Colorado, 80305, United States
## 73 253, Deer Trail Road, Lazy Acres, Boulder County, Colorado, 80455, United States
## 74 2430, Ginny Way, Blue Heron Estates, Lafayette, Boulder County, Colorado, 80026, United States
## 75 1313, Wildwood Court, Boulder, Boulder County, Colorado, 80305, United States
## 76 720, Yale Road, Boulder, Boulder County, Colorado, 80305, United States
## 78 7777, Country Creek Drive, Boulder, Niwot, Boulder County, Colorado, 80503, United States
## 80 3185, 29th Street, Boulder, Boulder County, Colorado, 80301, United States
## 81 770, South 45th Street, Boulder, Boulder County, Colorado, 80305, United States
## 82 6732, Snead Court, Boulder County, Colorado, 80503, United States
## 83 6430, Eagle Court, Boulder County, Colorado, 80503, United States
## 84 7157, Winona Court, Westminster, Adams County, Colorado, 80030, United States
## 85 1313, Wildwood Court, Boulder, Boulder County, Colorado, 80305, United States
## 86 Longmont, Boulder County, Colorado, 80503, United States
## 87 1052, Quince Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 88 1355, Redwood Avenue, Boulder, Boulder County, Colorado, 80304, United States
## 93 7971, Uptown Avenue, Arista, Broomfield, Colorado, 80021, United States
## 94 10081, Oak Street, Westminster, Jefferson County, Colorado, 80021, United States
## 95 4342, Augusta Drive, Broadlands, Broomfield, Colorado, 80023, United States
## 97 384 Driftwood Circle, 384, Driftwood Circle, Waneka Landing, Lafayette, Boulder County, Colorado, 80026, United States
## 99 Magic Mansion, 2267, Holyoke Drive, Boulder, Boulder County, Colorado, 80305, United States
## 101 3878, Hancock Drive, Boulder, Boulder County, Colorado, 80303, United States
## 103 Jay Road, Boulder, Boulder County, Colorado, 80301, United States
## 106 350, Bates Avenue, Boulder, Boulder County, Colorado, 80305, United States
## 115 850, Hartford Drive, Boulder, Boulder County, Colorado, 80305, United States
## 200 759, Nighthawk Circle, Louisville, Boulder County, Colorado, 80027, United States
## zip_code city_update
## 1 80304 Boulder
## 2 80503 Longmont
## 3 80305 Boulder
## 4 80503 Niwot
## 5 80305 Boulder
## 6 80504 Longmont
## 7 80305 Boulder
## 8 80303 Boulder
## 10 80301 Boulder
## 11 80301 Boulder
## 13 80544 Niwot
## 14 80032 States
## 16 80305 Boulder
## 17 80306 Boulder
## 18 80302 Boulder
## 19 80021 States
## 20 80304 Boulder
## 21 80305 Boulder
## 22 80027 Louisville
## 23 80020 States
## 24 80302 Boulder
## 25 80027 Louisville
## 26 80304 Lane
## 28 80304 Boulder
## 29 80305 Boulder
## 31 80303 Boulder
## 34 80027 Boulder
## 35 80304 Boulder
## 36 80027 Louisville
## 37 80305 Boulder
## 38 80304 Boulder
## 39 80305 Boulder
## 40 80504 Longmont
## 41 80304 Lane
## 42 80301 Boulder
## 43 80303 Boulder
## 44 80025 Springs
## 45 80027 Superior
## 46 80304 Hills
## 47 80027 Superior
## 48 80027 Boulder
## 49 80305 Boulder
## 50 80304 Boulder
## 51 80305 Boulder
## 52 80301 Boulder
## 53 80027 Louisville
## 55 80305 Boulder
## 56 2872 Lafayette
## 57 80027 Louisville
## 58 80027 Louisville
## 59 80304 Boulder
## 60 80303 Boulder
## 61 80301 Boulder
## 62 80302 Boulder
## 63 80020 States
## 64 80301 Boulder
## 65 80304 Hill
## 66 80026 Erie
## 67 80026 Lafayette
## 69 80027 Louisville
## 70 80302 Boulder
## 71 80301 Boulder
## 72 80305 Boulder
## 73 80455 Acres
## 74 80026 Lafayette
## 75 80305 Boulder
## 76 80305 Boulder
## 78 80503 Niwot
## 80 80301 Boulder
## 81 80305 Boulder
## 82 80503 Court
## 83 80503 Court
## 84 80030 States
## 85 80305 Boulder
## 86 80503 Longmont
## 87 80304 Boulder
## 88 80304 Boulder
## 93 80021 States
## 94 80021 States
## 95 80023 States
## 97 80026 Lafayette
## 99 80305 Boulder
## 101 80303 Boulder
## 103 80301 Boulder
## 106 80305 Boulder
## 115 80305 Boulder
## 200 80027 Louisville
## loc_status
## 1 outside - Complete data throughout fire period
## 2 outside - Complete data throughout fire period
## 3 outside - Complete data throughout fire period
## 4 outside - Complete data throughout fire period
## 5 outside - Sensor offline during fire, returned online
## 6 outside - Sensor offline during fire, did not return online
## 7 outside - Sensor offline during fire, returned online
## 8 outside - Sensor offline during fire, returned online
## 10 outside - Complete data throughout fire period
## 11 inside - Complete data throughout fire period
## 13 inside - Complete data throughout fire period
## 14 outside - Complete data throughout fire period
## 16 outside - Complete data throughout fire period
## 17 outside - Complete data throughout fire period
## 18 outside - Complete data throughout fire period
## 19 outside - Complete data throughout fire period
## 20 inside - Complete data throughout fire period
## 21 outside - Sensor offline during fire, returned online
## 22 outside - Sensor offline during fire, returned online
## 23 outside - Complete data throughout fire period
## 24 outside - Complete data throughout fire period
## 25 outside - Sensor offline during fire, returned online
## 26 outside - Sensor offline during fire, returned online
## 28 outside - Complete data throughout fire period
## 29 inside - Complete data throughout fire period
## 31 outside - Complete data throughout fire period
## 34 inside - Sensor offline during fire, returned online
## 35 inside - Sensor offline during fire, returned online
## 36 outside - Complete data throughout fire period
## 37 inside - Complete data throughout fire period
## 38 outside - Complete data throughout fire period
## 39 outside - Complete data throughout fire period
## 40 outside - Complete data throughout fire period
## 41 inside - Sensor offline during fire, returned online
## 42 inside - Complete data throughout fire period
## 43 inside - Complete data throughout fire period
## 44 outside - Sensor offline during fire, returned online
## 45 outside - Sensor offline during fire, did not return online
## 46 outside - Complete data throughout fire period
## 47 outside - Sensor offline during fire, returned online
## 48 outside - Sensor offline during fire, returned online
## 49 outside - Complete data throughout fire period
## 50 outside - Complete data throughout fire period
## 51 inside - Complete data throughout fire period
## 52 inside - Complete data throughout fire period
## 53 inside - Complete data throughout fire period
## 55 outside - Sensor offline during fire, returned online
## 56 outside - Complete data throughout fire period
## 57 outside - Complete data throughout fire period
## 58 outside - Complete data throughout fire period
## 59 outside - Complete data throughout fire period
## 60 outside - Complete data throughout fire period
## 61 outside - Complete data throughout fire period
## 62 outside - Complete data throughout fire period
## 63 outside - Complete data throughout fire period
## 64 inside - Complete data throughout fire period
## 65 inside - Sensor offline during fire, returned online
## 66 outside - Complete data throughout fire period
## 67 outside - Complete data throughout fire period
## 69 inside - Sensor offline during fire, returned online
## 70 inside - Sensor offline during fire, did not return online
## 71 outside - Complete data throughout fire period
## 72 inside - Sensor offline during fire, returned online
## 73 inside - Sensor offline during fire, returned online
## 74 inside - Complete data throughout fire period
## 75 inside - Complete data throughout fire period
## 76 inside - Sensor offline during fire, returned online
## 78 outside - Complete data throughout fire period
## 80 outside - Complete data throughout fire period
## 81 inside - Sensor offline during fire, returned online
## 82 outside - Sensor offline during fire, returned online
## 83 outside - Sensor offline during fire, returned online
## 84 outside - Complete data throughout fire period
## 85 outside - Complete data throughout fire period
## 86 outside - Complete data throughout fire period
## 87 inside - Complete data throughout fire period
## 88 outside - Complete data throughout fire period
## 93 outside - Complete data throughout fire period
## 94 inside - Complete data throughout fire period
## 95 outside - Complete data throughout fire period
## 97 outside - Complete data throughout fire period
## 99 outside - Sensor offline during fire, returned online
## 101 outside - Complete data throughout fire period
## 103 outside - Complete data throughout fire period
## 106 inside - Sensor offline during fire, returned online
## 115 inside - Sensor offline during fire, returned online
## 200 inside - Sensor offline during fire, did not return online# interactive map of damaged/destroyed buildings and AQ sensors (indoor/outdoor)
# chose not to do icons for destroyed/damaged buildings because the map got really busy
fire_destruction_and_AQ_plot = leaflet(sensors) %>%
addTiles() %>%
addAwesomeMarkers(icon = ~sensorIcon[Location],
lng = sensors$Lon, lat = sensors$Lat) %>%
addCircleMarkers(color = ~pal(destroyed_damaged$type),
radius = 3.5,
opacity = 1,
lng = destroyed_damaged$long, lat = destroyed_damaged$lat) %>%
addLegend("topright", pal = pal, values = ~destroyed_damaged$type,
title = "Fire damage")Leaflet map that colors each sensor icon based on the time period it has data:
# create palette for location + time period
timeSensorIcon = awesomeIconList(
"inside - Complete data throughout fire period" = makeAwesomeIcon(
icon = "home",
iconColor = "white",
markerColor = "green",
library = "fa"
),
"inside - Sensor offline during fire, did not return online" = makeAwesomeIcon(
icon = "tree",
iconColor = "white",
markerColor = "gray",
library = "fa"
),
"inside - Sensor offline during fire, returned online" = makeAwesomeIcon(
icon = "home",
iconColor = "white",
markerColor = "blue",
library = "fa"
),
"outside - Complete data throughout fire period" = makeAwesomeIcon(
icon = "home",
iconColor = "white",
markerColor = "green",
library = "fa"
),
"outside - Sensor offline during fire, did not return online" = makeAwesomeIcon(
icon = "tree",
iconColor = "white",
markerColor = "gray",
library = "fa"
),
"outside - Sensor offline during fire, returned online" = makeAwesomeIcon(
icon = "tree",
iconColor = "white",
markerColor = "blue",
library = "fa"
))
# map based on when sensor was added
(fire_destruction_and_AQ_plot = leaflet(sensors) %>%
addTiles() %>%
addAwesomeMarkers(icon = ~timeSensorIcon[loc_status],
lng = sensors$Lon, lat = sensors$Lat) %>%
addCircleMarkers(color = ~pal(destroyed_damaged$type),
radius = 3.5,
opacity = 1,
lng = destroyed_damaged$long, lat = destroyed_damaged$lat) %>%
addLegend("topright", pal = pal, values = ~destroyed_damaged$type,
title = "Fire damage") %>%
addLegend("bottomright",
pal = colorFactor(c("green","gray","blue"), domain = unique(sensors$Status)),
values = ~sensors$Status, title = "Sensor Time Info"))Create spatial objects
To do more in-depth spatial analysis and mapping, we need to convert our data into spacetime objects.
## Create spacetime objects -- 10 minute averages
# first, re-format datetime as a xts object
AQ_df$datetime = as.POSIXct(AQ_df$datetime)
# create spatial points object for each sensor
PA_sensors = SpatialPoints(AQ_df[!duplicated(AQ_df$ID), c("Lon", "Lat")],proj4string = CRS(prg))
summary(PA_sensors)## Object of class SpatialPoints
## Coordinates:
## min max
## Lon -105.35808 -105.02234
## Lat 39.82646 40.13818
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Number of points: 144# construct spatiotemporal object
PA_STFDF = stConstruct(AQ_df, space = c("Lon","Lat"), time = "datetime", crs = CRS(prg), SpatialObj = sensors)
# turn ST object into a STFDF (stationary points)
PA_STFDF = as(PA_STFDF,"STFDF")
# see summary of data
summary(PA_STFDF)## Object of class STFDF
## with Dimensions (s, t, attr): (144, 17706, 12)
## [[Spatial:]]
## Object of class SpatialPoints
## Coordinates:
## min max
## Lon -105.35808 -105.02234
## Lat 39.82646 40.13818
## Is projected: FALSE
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Number of points: 144
## [[Temporal:]]
## Index timeIndex
## Min. :2021-12-30 00:00:00 Min. : 1
## 1st Qu.:2022-01-29 17:42:30 1st Qu.: 4427
## Median :2022-03-01 11:25:00 Median : 8854
## Mean :2022-03-01 11:25:00 Mean : 8854
## 3rd Qu.:2022-04-01 06:07:30 3rd Qu.:13280
## Max. :2022-05-01 23:50:00 Max. :17706
## [[Data attributes:]]
## X ID Name Location
## Min. : 1 Min. : 3062 Length:2549664 Length:2549664
## 1st Qu.: 482877 1st Qu.: 82879 Class :character Class :character
## Median : 965752 Median :110042 Mode :character Mode :character
## Mean : 965752 Mean : 99993
## 3rd Qu.:1448628 3rd Qu.:128801
## Max. :1931503 Max. :146952
## NA's :618161 NA's :618161
## pm25_a pm25_b temp rh
## Min. : 0.0 Min. : 0.0 Min. :-234.0 Min. : 0.0
## 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 38.2 1st Qu.:19.0
## Median : 0.5 Median : 0.4 Median : 50.0 Median :27.4
## Mean : 34.6 Mean : 6.6 Mean : 50.2 Mean :30.0
## 3rd Qu.: 3.9 3rd Qu.: 3.7 3rd Qu.: 62.6 3rd Qu.:40.4
## Max. :4999.0 Max. :3542.5 Max. : 127.4 Max. :94.6
## NA's :618161 NA's :618161 NA's :652435 NA's :652435
## a_b_agree pm_avg ab_difference corrected_pm
## Length:2549664 Min. : 0.0 Min. :-3537.2 Min. : -2.0
## Class :character 1st Qu.: 0.0 1st Qu.: -0.1 1st Qu.: 3.3
## Mode :character Median : 0.6 Median : 0.0 Median : 4.2
## Mean : 20.6 Mean : 28.0 Mean : 18.6
## 3rd Qu.: 4.0 3rd Qu.: 0.2 3rd Qu.: 5.4
## Max. :2514.5 Max. : 4994.6 Max. :1978.1
## NA's :618161 NA's :618161 NA's :652435# object class
class(PA_STFDF)## [1] "STFDF"
## attr(,"package")
## [1] "spacetime"# object dimensions
dim(PA_STFDF)## space time variables
## 144 17706 12Map of all sensors we have data for (many not online during the fire)
ggplot(all_bounds, aes(geometry=geometry)) +
geom_sf() +
geom_sf(data=marshall_fire$geometry, fill="red", alpha=0.6) +
geom_sf(data=st_as_sf(PA_sensors))
Time series plots
## Time series plots
stplot(PA_STFDF[,"2021-12-30::2021-12-31","temp"],mode="ts")
stplot(PA_STFDF[,"2021-12-30::2021-12-31","rh"],mode="ts")
stplot(PA_STFDF[,"2021-12-30::2021-12-31","corrected_pm"],mode="ts")
Spacetime Plots
These show time on the y-axis, each sensor on the x-axis, and the color corresponds to the variable (temperature, relative humidity, or PM2.5).
stplot(PA_STFDF[,"2021-12-30::2021-12-31","temp"],mode="xt")
stplot(PA_STFDF[,"2021-12-30::2021-12-31","rh"],mode="xt")
stplot(PA_STFDF[,"2021-12-30::2021-12-31","corrected_pm"],mode="xt")
Sensors With Complete Data
# make an STFDF with only the sensors that had complete data for the fire
complete_sensors = sensors %>% filter(Status == "Complete data throughout fire period")
complete_STFDF = PA_STFDF[PA_STFDF@data$ID %in% complete_sensors$ID]
stplot(complete_STFDF[,,"corrected_pm"], mode="ts")
## Plotting the PM2.5 values for each sensor on 12/30/2021
stplot(complete_STFDF[,"2021-12-30::2021-12-31","corrected_pm"], mode="ts")
stplot(complete_STFDF[,"2021-12-30::2021-12-31","corrected_pm"], mode="tp")
Kriging
Funky Data from Variograms
After running through this kriging methodology without removing any data points, there were a few sensors that had very strange values when looking through the variograms. This prompted me to plot them individually and examine the actual PM2.5 values for the fire period, to see if any of them seemed unreasonable. Both the Eisenhower Dr and 1333 Wildwood Ct sensors have values that are far too high, so those are dropped for our analysis done below.
# Index 15 -> ID 60517, STANDLEY LAKE
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 60517]["2021-12-30::2022-01-02", "corrected_pm"]))
# looks normal
# Index 23 -> ID 91877, Eisenhower Dr
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 91877]["2021-12-30::2022-01-02", "corrected_pm"]))
# has incredibly large values midday Friday - Sunday, could be snow related?
# Index 36 -> ID 112606, 1333 Wildwood Ct., Boulder, CO 80305 sensor
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 112606]["2021-12-30::2022-01-02", "corrected_pm"]))
# looking at purpleair map, A channel is 'downgraded', 0% confidence
# Index 40 -> ID 119547, Broadlands sensor
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 119547]["2021-12-30::2022-01-02", "corrected_pm"]))
# reasonable valuesLooking at Dec 30 - Jan 2, only outside sensors.
# get all of the data for the 12/30 to 1/2 timeframe
krig_data <- complete_STFDF[,"2021-12-30::2022-01-02"]
# limit to only outdoor sensors
krig_data <- krig_data[krig_data@data$Location == "outside"]
# drop sensor 112606 (Wildwood Ct) because it has 0% confidence according to PurpleAir map
krig_data <- krig_data[krig_data@data$ID != 112606]
# drop values for "Eisenhower Dr" sensor that are over 1000
krig_data@data[(krig_data@data$ID == 91877) &
(as.numeric(krig_data@data$corrected_pm) > 1000) &
!is.na(krig_data@data$corrected_pm),]$corrected_pm <- NA
# get rid of the weird sensor that might be throwing residuals off
krig_data <- krig_data[krig_data@data$Name != "ntsky"]
# reproject the data to utm
reproj <- "+proj=utm +zone=13 +datum=WGS84 +units=km"
krig_data@sp <- spTransform(krig_data@sp, reproj)
# aggregate from 10 min intervals to 1 hour intervals
krig_data <- aggregate(krig_data, by="hour", mean)
# split the data into the days we want for the fire period
dec30 <- krig_data[, "2021-12-30"]
dec31 <- krig_data[, "2021-12-31"]
jan1 <- krig_data[, "2022-01-01"]
jan2 <- krig_data[, "2022-01-02"]
# split each day into the time periods of focus
# pre-fire period: betweeen 00:00 - 11:00 MST on 12/30/2021
pre_fire <- aggregate(dec30, by="11 hours", mean)[, "2021-12-30 00:00:00"]
# fire period (11am - 5pm)
during_fire <- aggregate(dec30[, "2021-12-30 11:00:00::2021-12-30 16:00:00"], by="6 hours", mean)
# post-fire on 12/30/2021
evening_of_fire <- aggregate(dec30[, "2021-12-30 16:00:00::2021-12-31 00:00:00"], by="8 hours", mean)
# dec 31st
december31 <- aggregate(dec31, by="day", mean)
# january 1st
january1 <- aggregate(jan1, by="day", mean)
# january 2nd
january2 <- aggregate(jan2, by="day", mean)
periods <- c(pre_fire, during_fire, evening_of_fire, december31, january1, january2)
# create variograms for these time periods
pre_fire.vgm <- variogram(corrected_pm~1, pre_fire[!is.na(pre_fire$corrected_pm),])
during_fire.vgm <- variogram(corrected_pm~1, during_fire[!is.na(during_fire$corrected_pm),])
evening_of_fire.vgm <- variogram(corrected_pm~1, evening_of_fire[!is.na(evening_of_fire$corrected_pm),])
december31.vgm <- variogram(corrected_pm~1, december31[!is.na(december31$corrected_pm),])
january1.vgm <- variogram(corrected_pm~1, january1[!is.na(january1$corrected_pm),])
january2.vgm <- variogram(corrected_pm~1, january2[!is.na(january2$corrected_pm),])
# plot the variograms
plot(pre_fire.vgm)
plot(during_fire.vgm)
plot(evening_of_fire.vgm)
plot(december31.vgm)
plot(january1.vgm)
plot(january2.vgm)
Variogram Clouds
for (period in periods) {
cloud <- variogram(corrected_pm~1, period[!is.na(period$corrected_pm), ], cloud=TRUE)
map <- variogram(corrected_pm~1, period[!is.na(period$corrected_pm), ], map=TRUE, cutoff=20, width=2)
print(plot(cloud))
print(plot(map))
}
No good variogram models!!!!
We’re going to try autofitting variograms & models to the variograms to see if this returns anything different than above.
pre_fire.fitted <- fit.variogram(pre_fire.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fitpre_fire.fitted## model psill range
## 1 Nug 86.24564 0.0
## 2 Exp 62387.14334 -10613.7plot(pre_fire.vgm, model=pre_fire.fitted)
## model psill range
## 1 Nug 86.24564 0.0
## 2 Exp 62387.14334 -10613.7during_fire.fitted <- fit.variogram(during_fire.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?during_fire.fitted## model psill range
## 1 Nug 152.5208 0.00000
## 2 Exp 938.8172 68.95887plot(during_fire.vgm, model=during_fire.fitted)
evening_of_fire.fitted <- fit.variogram(evening_of_fire.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?evening_of_fire.fitted## model psill range
## 1 Nug 0.0000 0.00000
## 2 Gau 228.5243 23.99524plot(evening_of_fire.vgm, model=evening_of_fire.fitted)
december31.fitted <- fit.variogram(december31.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?december31.fitted## model psill range
## 1 Nug 1.141338 0.00000
## 2 Gau 38.915184 16.75872plot(december31.vgm, model=december31.fitted)
january1.fitted <- fit.variogram(january1.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?january1.fitted## model psill range
## 1 Nug 4.897684 0.000000
## 2 Sph 7.572974 2.524618plot(january1.vgm, model=january1.fitted)
january2.fitted <- fit.variogram(january2.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
january2.fitted## model psill range
## 1 Nug 1.572919 0.00000
## 2 Exp 33.290787 12.57439plot(january1.vgm, model=january2.fitted)
These models are no good (some can’t even find models), so kriging isn’t going to be working with this data.
Breaking down into smaller timeframes
Let’s see if things look different when we examine each hour, instead of aggregating into ~6 hour periods.
for (i in 0:11) {
if (i < 10) {
i <- paste("0", i, sep="")
}
t <- paste("2021-12-30 ", i, ":00:00", sep="")
i.data <- dec30[, t]
i.vgm <- variogram(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),])
# print(plot(i.vgm, main=paste("Variogram at", t)))
i.fitted <- fit.variogram(i.vgm, vgm(c("Exp", "Sph", "Gau", "Mat")))
print(i.fitted)
print(plot(i.vgm, model=i.fitted, main=paste("Variogram at", t)))
}## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## model psill range
## 1 Nug 0.000000 0.000000
## 2 Sph 3.904079 5.151113## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.000000 0.000000
## 2 Sph 3.095985 3.763504## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.000000 0.000000
## 2 Gau 3.260366 1.309588## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.4409823 0.000000
## 2 Sph 1.0171682 3.240835## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.1678131 0.00000
## 2 Gau 0.2917075 10.77585## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.0501557 0.00000
## 2 Gau 1.7184797 71.28375## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.028565 0.00000
## 2 Gau 0.201726 17.88329## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : linear model has singular covariance matrix
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : linear model has singular covariance matrix## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.01918062 0.00000
## 2 Sph 0.01776751 1.34037## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 0.034425382 0.00000
## 2 Sph 0.005219032 29.42497## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## model psill range
## 1 Nug 6181.04 0.0000
## 2 Exp 50155.36 -305.7408## model psill range
## 1 Nug 6181.04 0.0000
## 2 Exp 50155.36 -305.7408## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit
## model psill range
## 1 Nug 391.9872 0.0000
## 2 Exp 375.1265 -189.0302## model psill range
## 1 Nug 391.9872 0.0000
## 2 Exp 375.1265 -189.0302## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, model, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : singular model in variogram fit## Warning in fit.variogram(object, x, fit.sills = fit.sills, fit.ranges =
## fit.ranges, : No convergence after 200 iterations: try different initial values?
## model psill range
## 1 Nug 1.259568 0.000000
## 2 Gau 0.000000 4.613097
Still no good variograms!
Making grid for interpolation
We create a 50x50 grid for interpolation. This results in the resolution being 0.009372036 km wide by 0.007888942 km tall.
grd <- SpatialPixels(SpatialPoints(as_Spatial(st_make_grid(all_bounds, n=50))), proj4string = prg)
plot(grd)
grd <- grd[as_Spatial(all_bounds),]
plot(grd)
grd <- spTransform(grd, CRS(reproj))Plotting Corrected PM
This is just used to check the data from the sensors, and visualize if any of the sensors have odd values.
plot(pre_fire$corrected_pm)
plot(during_fire$corrected_pm)
plot(evening_of_fire$corrected_pm)
plot(december31$corrected_pm)
plot(january1$corrected_pm)
plot(january2$corrected_pm)
IDW
Inverse distance weighting interpolation is used to get a map of values across our new grid object for each time period.
pre_fire.idw <- idw(corrected_pm~1, pre_fire[!is.na(pre_fire$corrected_pm),], grd)## [inverse distance weighted interpolation]spplot(pre_fire.idw[,1], main="Pre-Fire Period on 12/30")
during_fire.idw <- idw(corrected_pm~1, during_fire[!is.na(during_fire$corrected_pm),], grd)## [inverse distance weighted interpolation]spplot(during_fire.idw[,1], main="During Fire on 12/30")
evening_of_fire.idw <- idw(corrected_pm~1, evening_of_fire[!is.na(evening_of_fire$corrected_pm),], grd)## [inverse distance weighted interpolation]spplot(evening_of_fire.idw[,1], main="After Fire on 12/30")
december31.idw <- idw(corrected_pm~1, december31[!is.na(december31$corrected_pm),], grd)## [inverse distance weighted interpolation]spplot(december31.idw[,1], main="December 31st")
january1.idw <- idw(corrected_pm~1, january1[!is.na(january1$corrected_pm),], grd)## [inverse distance weighted interpolation]spplot(january1.idw[,1], main="January 1st")
january2.idw <- idw(corrected_pm~1, january2[!is.na(january2$corrected_pm),], grd)## [inverse distance weighted interpolation]spplot(january2.idw[,1], main="January 2nd")
Let’s make the maps easier to read. We use the EPA AQI mapping standards for these and create plots that have consistent legends.
aqi_colors <- c("#00e400",
"#ffff00",
"#ff7e00",
"#ff0000",
"#8f3f97",
"#7e0023")
aqi_labels <- c("Good (0-12.0)",
"Moderate (12.1-35.4)",
"Unhealthy for Sensitive Groups (35.5-55.4)",
"Unhealthy (55.5-150.4)",
"Very Unhealthy (150.5-250.4)",
"Hazardous (250.5-500.4)")
factorize_aqi <- function (idw) {
idw$PM2.5 <- cut(idw$var1.pred, breaks=c(0, 12.0, 35.4, 55.4, 150.4, 250.4), na.omit=T)
return(idw)
}
plot_idw <- function(idw, title) {
return(ggplot(all_bounds) +
geom_sf() +
geom_sf(data=st_as_sf(factorize_aqi(idw)), aes(col=PM2.5), size=2) +
scale_color_manual(values=aqi_colors, labels=aqi_labels, name="PM2.5 (µg/m^3)") +
geom_sf(data=marshall_fire$geometry, col="gray20", alpha=0.1) +
ggtitle(title))
}
plot_idw(pre_fire.idw, "Pre Fire IDW")
plot_idw(during_fire.idw, "During Fire IDW")
plot_idw(evening_of_fire.idw, "Evening of fire IDW")
plot_idw(december31.idw, "December 31st IDW")
plot_idw(january1.idw, "January 1st IDW")
plot_idw(january2.idw, "January 2nd IDW")
Now we take the above plots & plot them as rasters to clear up any visual confusion.
r_grd <- raster(grd, nrows=50, ncols=49) # creating a 50x50 object gives some weird visuals
plot_idw_raster <- function(idw, title) {
r <- raster::rasterize(idw, r_grd, field="var1.pred")
extent(r) <- st_bbox(all_bounds)
r_df <- r %>%
rasterToPoints() %>%
data.frame() %>%
mutate(cuts = cut(layer, breaks=c(0, 12.0, 35.4, 55.4, 150.4, 250.4), na.omit=T))
return(ggplot() +
geom_tile(data=r_df, aes(x=x, y=y, fill=cuts))+
scale_fill_manual(values=aqi_colors, labels=aqi_labels, name="PM2.5 (µg/m^3)", drop=FALSE) +
geom_sf(data=marshall_fire$geometry, col="gray20", alpha=0.1) +
ggtitle(title))
}
plot_idw_raster(pre_fire.idw, "Pre Fire IDW")
plot_idw_raster(during_fire.idw, "During Fire IDW")
plot_idw_raster(evening_of_fire.idw, "Evening of fire IDW")
plot_idw_raster(december31.idw, "December 31st IDW")
plot_idw_raster(january1.idw, "January 1st IDW")
plot_idw_raster(january2.idw, "January 2nd IDW")
Export smoke affected area
r <- raster::rasterize(during_fire.idw, r_grd, field="var1.pred")
extent(r) <- st_bbox(all_bounds)
r[r < 12] <- NA
writeRaster(r, "../GIS_inputs_destruction_fireboundary/smoke_affected.tiff", overwrite=TRUE)Cross Validation
Using leave-one-out cross validation to get residuals for the sensors for each time period. Each sensor is left out once, then the average residual is calculated from the other 37 times the model is created.
pre_fire.idw_cv <- krige.cv(corrected_pm~1, pre_fire[!is.na(pre_fire$corrected_pm),], nfold=38)
pre_fire.idw_sd <- sd(pre_fire.idw_cv$residual)
paste("Pre fire standard deviation:", pre_fire.idw_sd)## [1] "Pre fire standard deviation: 7.58649993004707"spplot(pre_fire.idw_cv, main="Pre-Fire Period on 12/30")
pre_fire.idw_cv$zscore <- pre_fire.idw_cv$residual / pre_fire.idw_sd
during_fire.idw_cv <- krige.cv(corrected_pm~1, during_fire[!is.na(during_fire$corrected_pm),], nfold=38)
during_fire.idw_sd <- sd(during_fire.idw_cv$residual)
paste("During fire standard deviation:", during_fire.idw_sd)## [1] "During fire standard deviation: 15.376531044756"spplot(during_fire.idw_cv, main="During Fire on 12/30")
during_fire.idw_cv$zscore <- during_fire.idw_cv$residual / during_fire.idw_sd
evening_of_fire.idw_cv <- krige.cv(corrected_pm~1, evening_of_fire[!is.na(evening_of_fire$corrected_pm),], nfold=38)
evening_of_fire.idw_sd <- sd(evening_of_fire.idw_cv$residual)
paste("Evening of fire standard deviation:", evening_of_fire.idw_sd)## [1] "Evening of fire standard deviation: 4.85088863186492"spplot(evening_of_fire.idw_cv, main="After Fire on 12/30")
evening_of_fire.idw_cv$zscore <- evening_of_fire.idw_cv$residual / evening_of_fire.idw_sd
december31.idw_cv <- krige.cv(corrected_pm~1, december31[!is.na(december31$corrected_pm),], nfold=38)
december31.idw_sd <- sd(december31.idw_cv$residual)
paste("December 31st standard deviation:", december31.idw_sd)## [1] "December 31st standard deviation: 2.36672112753216"spplot(december31.idw_cv, main="December 31st")
december31.idw_cv$zscore <- december31.idw_cv$residual / december31.idw_sd
january1.idw_cv <- krige.cv(corrected_pm~1, january1[!is.na(january1$corrected_pm),], nfold=38)
january1.idw_sd <- sd(january1.idw_cv$residual)
paste("January 1 standard deviation:", january1.idw_sd)## [1] "January 1 standard deviation: 3.68472999489718"spplot(january1.idw_cv, main="January 1st")
january1.idw_cv$zscore <- january1.idw_cv$residual / january1.idw_sd
january2.idw_cv <- krige.cv(corrected_pm~1, january2[!is.na(january2$corrected_pm),], nfold=38)
january2.idw_sd <- sd(january2.idw_cv$residual)
paste("January 2 standard deviation:", january2.idw_sd)## [1] "January 2 standard deviation: 3.6116917414784"spplot(january1.idw_cv, main="January 2nd")
january2.idw_cv$zscore <- january2.idw_cv$residual / january2.idw_sdPlotting these residuals & normalizing them to z-scores for easier visual comparison.
plot_idw_resid <- function(idw.cv, title) {
return(ggplot(all_bounds) +
geom_sf(fill="#FbFbFb", color="gray") +
geom_sf(data=marshall_fire$geometry, col="red", alpha=0.1) +
geom_sf(data=st_as_sf(idw.cv), aes(fill=residual), color="black", shape=21, size=2) +
scale_color_gradient2(aesthetics="fill") +
ggtitle(title))
}
plot_idw_resid_zscore <- function(idw.cv, title) {
return(ggplot(all_bounds) +
geom_sf(fill="#FbFbFb", color="gray") +
geom_sf(data=marshall_fire$geometry, col="red", alpha=0.1) +
geom_sf(data=st_as_sf(idw.cv), aes(fill=zscore), color="black", shape=21, size=2) +
scale_color_gradient2(aesthetics="fill") +
ggtitle(title))
}
plot_idw_resid(pre_fire.idw_cv, "Pre-Fire Residuals")
plot_idw_resid(during_fire.idw_cv, "During Fire Residuals")
plot_idw_resid(evening_of_fire.idw_cv, "After Fire Residuals")
plot_idw_resid(december31.idw_cv, "December 31st Residuals")
plot_idw_resid(january1.idw_cv, "January 1st Residuals")
plot_idw_resid(january2.idw_cv, "January 2nd Residuals")
plot_idw_resid_zscore(pre_fire.idw_cv, "Pre-Fire Residual Z-Scores")
plot_idw_resid_zscore(during_fire.idw_cv, "During Fire Residual Z-Scores")
plot_idw_resid_zscore(evening_of_fire.idw_cv, "After Fire Residual Z-Scores")
plot_idw_resid_zscore(december31.idw_cv, "December 31st Residual Z-Scores")
plot_idw_resid_zscore(january1.idw_cv, "January 1st Residual Z-Scores")
plot_idw_resid_zscore(january2.idw_cv, "January 2nd Residual Z-Scores")
Trying to plot each map with the residuals next to it, however these just look messy.
plot_idw_cv <- function(idw, idw_title, idw.cv) {
p1 <- plot_idw(idw, idw_title)
p2 <- plot_idw_resid(idw.cv, "Residuals")
return(grid.arrange(p1, p2, ncol=2))
}
plot_idw_cv(pre_fire.idw, "Pre-Fire IDW", pre_fire.idw_cv)
plot_idw_cv(during_fire.idw, "During Fire IDW", during_fire.idw_cv)
plot_idw_cv(evening_of_fire.idw, "Evening of Fire IDW", evening_of_fire.idw_cv)
plot_idw_cv(december31.idw, "December 31st IDW", december31.idw_cv)
plot_idw_cv(january1.idw, "January 1st IDW", january1.idw_cv)
plot_idw_cv(january2.idw, "January 2nd IDW", january2.idw_cv)
Notes:
- Pre-fire period the residuals look fairly good, some over & underprediction in downtown Boulder
- During fire period has 1 large underprediction value, few other overpredicted near fire area
- Evening after the fire has underprediction, but not much
- December 31st has high under & over predicted values right next to each other near fire
- Same as December 31st on Jan 1st
- Same as prior 2 days
- Something happened starting on December 31st with these two sensors? Will look into both and see how the values differ over time
# One sensor is "Eisenhower Dr"
plot(fortify(PA_STFDF[PA_STFDF@data$Name == "Eisenhower Dr"]["2021-12-30::2022-01-02", "corrected_pm"]))
# The other is "Purple House"
plot(fortify(PA_STFDF[PA_STFDF@data$Name == "Purple House"]["2021-12-30::2022-01-02", "corrected_pm"]), col="red")
The issue is coming from the Eisenhower Dr sensor most likely, as this is the one that seems messed up by the snow. Going to see how everything looks excluding the high values for this sensor. This fed back into our analysis, and has already been accounted for in the current models.
plot(fortify(PA_STFDF[PA_STFDF@data$Name == "Eisenhower Dr"]["2021-12-30::2022-01-02", "corrected_pm"]) %>%
filter(corrected_pm <= 1000))
Checking the sensor on the edge of the fire boundary with high residuals/z scores.
plot(fortify(PA_STFDF[PA_STFDF@data$Name == "ntsky"]["2021-12-30::2022-01-02", "corrected_pm"]))
Looping at Hourly Timescale
We create gifs to make visualization easier, as well as looking at hourly timescales instead of the six predefined fire periods.
for (j in c(30, 31, 01, 02)) {
if (j > 3) {
date <- paste("2021-12-", j, sep="")
} else {
date <- paste("2022-01-", j, sep="")
}
for (i in 0:23) {
if (i < 10) {
i <- paste("0", i, sep="")
}
t <- paste(date, " ", i, ":00:00", sep="")
i.data <- krig_data[, t]
i.idw <- idw(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),], grd)
print(plot_idw_raster(i.idw, paste(t, "IDW")))
}
}
for (j in c(30, 31, 01, 02)) {
if (j > 3) {
date <- paste("2021-12-", j, sep="")
} else {
date <- paste("2022-01-", j, sep="")
}
for (i in 0:23) {
if (i < 10) {
i <- paste("0", i, sep="")
}
t <- paste(date, " ", i, ":00:00", sep="")
i.data <- krig_data[, t]
i.idw <- idw(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),], grd)
i.idw_cv <- krige.cv(corrected_pm~1, i.data[!is.na(i.data$corrected_pm),], nfold=38)
i.idw_sd <- sd(i.idw_cv$residual)
i.idw_cv$zscore <- i.idw_cv$residual / i.idw_sd
print(plot_idw_resid_zscore(i.idw_cv, paste(t, "IDW Residual Z-Scores")))
}
}
Daily Gifs Starting 1/3/2022
# get the full timerange of data & clean it up
daily <- complete_STFDF[,"2022-01-01::"]
daily <- daily[daily@data$Location == "outside"]
# drop sensor 112606 b/c 0% confidence according to PurpleAir
daily <- daily[daily@data$ID != 112606]
# drop values for "Eisenhower Dr" sensor that're over 1000
daily@data[(daily@data$ID == 91877) &
(as.numeric(daily@data$corrected_pm) > 1000) &
!is.na(daily@data$corrected_pm),]$corrected_pm <- NA
# get rid of the weird sensor that might be throwing residuals off
daily <- daily[daily@data$Name != "ntsky"]
# get rid of "ponderosa" sensor values after april 11 (when sensor looks to get wonky)
daily@data[(daily@data$ID == 102614) & (daily@endTime < "2022-04-11") & !is.na(daily@data)] <- NA
# reproject & aggregate
daily@sp <- spTransform(daily@sp, reproj)
daily <- aggregate(daily, by="day", mean)for (j in 1:4) {
date <- paste("2022-0", j, sep="")
if (j %in% c(1, 3, 5)) {
range = 1:31
} else if (j == 2) {
range = 1:28
} else {
range = 1:30
}
for (i in range) {
if (i < 10) {
i <- paste("0", i, sep="")
}
day <- paste(date, "-", i, sep="")
data <- daily[, day]
idw <- idw(corrected_pm~1, data[!is.na(data$corrected_pm),], grd)
print(plot_idw_raster(idw, paste(day, "IDW")))
}
}
Weird Sensor from Gif
This sensor seemed to be really skewing our results starting in April 2022. Looking at the purple air website, this sensor has 0% confidence now, with values looking to lose credibility around the start of April, so earlier all of those values were dropped. (this already fed back into the current analysis we’re looking at).
plot(daily[, "2022-04-25"]$corrected_pm) # sensor at index 29 is weird = sensor ID 102614
plot(fortify(PA_STFDF[PA_STFDF@data$ID == 102614]['2022-04', "corrected_pm"]))
# ponderosa sensor has 0% confidence on purpleair map at present date